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Introduction 


This final report summarizes three years of work characterizing the radiometry of the 
Landsat 4, 5 and 7 Thematic Mappers. It is divided into six sections that are 
representative of the major areas of effort: 

Internal Calibrator Lamp Monitoring 
Vicarious Calibration 
Relative Gain Analysis 
Outgassing 

Landsat 4 Absolute Calibration 
Landsat 5 Scene Invariant Analysis. 

Each section provides a summary overview of the work that has been performed at 
SDSU. Major results are highlighted. In several cases, references are given to 
publications that have developed from this work. Several team members contributed to 
this report: Tim Ruggles, Dave Aaron, Shriharsha Madhavan, Esad Micijevic, Cory 
Mettler, and Jim Dewald. 

At the end of the report is a summary section. 
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LANDSAT 5 TM 1C RESPONSE 


Introduction 

This section of the final report summarizes the trends observed in the Landsat 5 
IC response data. It is currently believed that the changes occurring later in the 
instrument lifetime, and shown in the following plots, are most likely due to significant 
changes in one or more components within the IC subsystem, rather than changes to the 
detector arrays that would affect their response. This hypothesis tends to be supported by 
the fact that overall image quality, even in bumper mode operation, remains very high. 


Review of Data Processing Methodology 

The data under consideration here consist of Level 0 calibration data 
systematically corrected for scan-correlated-shift (SCS) artifacts. Data from Bands 1-4 
were also corrected for memory effect (ME) artifacts. Bias-corrected or ‘net’ integrated 
pulse values (NPV’s) were then extracted for each detector from the artifact-corrected IC 
data. After extraction, the NPV’s were further corrected for temperature effects (Bands 
1-4) and outgassing effects (Bands 5 and 7). Finally, a band-level NPV was estimated by 
averaging the corrected NPV’s across all operative detectors (due to an apparent error in 
the TMACS system that treated detector 4 in Band 2 and detector 3 in Band 5 as ‘dead’ — 
as is the case for the Landsat 4 TM — NPV’s from these detectors were not included in the 
calculations). 


Lifetime IC Response in LS [001] and LS [100], Bands 1-4 

Prior to year 1990, abrupt changes in IC response, on the order of 5% to 10% 
depending on the band, have been observed in Bands 1-4 of the individual lamp states 
[001] and [100]. Figures la and lb show these changes observed in the reverse scan data 
for Band 1 . Due to these early discontinuities, they were not seriously considered for 
future IC-based calibration efforts. 
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Figure la. Band 1 [001 1 IC Lifetime Response. 
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Figure lb. Band 1 [100] Lifetime IC Response. 


Several significant response ‘spikes’ were also observed for these lamp states, 
beginning in 1998 for LS [100] and 2000 for LS [001], indicating the occurrence of 
intermittent ‘flare-ups’ in the lamps. These changes have continued to occur after the 
transition to bumper mode operation in early 2002, as shown in the Band 1 plots of 
Figures 2a and 2b. For LS [001], the maximum variation observed in this period is on the 
order of 10% to 16%, depending on the band. There are periods of relatively stable 
response prior to and following these spikes. There also appears to be a general decrease 
in response, on the order of 10% to 16%, observed in LS [100] during this same period. 
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Figure 2a. Band 1 |001] Bumper Mode IC Response. Two spikes in response occurred in 2003; the 
larger spike in mid-2003 was on the order of 10-15% (depending on band). Overall lamp output 
appeared to be relatively stable prior to and following these spikes. Smaller response spikes may 
have occurred in early 2004, followed by even smaller spikes later on. 
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Figure 2b. Band 1 [100] Bumper Mode IC Response. An abrupt change in response was observed in 
mid 2002, following the transition to bumper mode operation. In addition to apparent response 
spikes occurring during early to mid 2003, there appeared to be a general decrease in response on the 
order of 10% to 15% during this period. Note the periods of relatively stable response from mid 
2002 to early 2003 and mid 2003 after each response spike; another period of stable response was 
observed after a slight response dip in mid-2004. 
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Lifetime 1C Response in LS[001] and LS[100], Bands 5 and 7 

Figures 3a and 3 b show the lifetime IC response in Band 5 for LS [001] and LS 
[100], respectively. 



Figure 3a. Band 5 [001] Lifetime IC Response. 



Figure 3b. Band 5 [100| Lifetime IC Response. 


Several things are noticeable from these plots. First, the discontinuities in Band 
1-4 responses observed in 1986 and 1988 cannot be readily observed in these bands. 
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Lifetime 1C Response in LS [010], Bands 1-4 

The remaining lamp state, [010], exhibited no apparent discontinuities in response 
in any reflective band prior to year 2000, as can be seen from the Band 1 plot shown in 
Figure 5. A reasonably simple ‘normalized’ model was derived from these data sets, 
which consisted of an early exponential decrease in response (believed to be due to 
outgassing from the spectral filters) followed by a linear increase in response. In general, 
these models were able to account for up to 95% of the observed variability in response; 
however, a second-order long period oscillatory effect could not be accounted for. The 
later linear increase was judged to be a ‘false’ effect, more likely due to changes in 
various components within the IC subsystem rather than changes to the detectors 
themselves; this hypothesis was supported by the fact that no such response trend had 
been observed in the vicarious calibration data, which are acquired through a different 
optical pathway within the instrument. As a result, the linear term from the normalized 
models was removed, leading to a normalized ‘exponential-only’ model to characterize 
the IC response. A lifetime calibration model for each band was obtained by scaling the 
exponential-only models to the estimated gain obtained in a cross-calibration 
measurement with the Landsat 7 ETM+ performed June 1-4, 1999. The models are 
currently implemented in a lookup table used by the NLAPS processing system. 



Figure 5. Band 1 (010] Lifetime IC Response. 


Following an apparent increase in the response beginning in 2000, a significant 
response spike was observed to occur in mid to late 2001; this post-2000 behavior 
effectively rendered extension of the calibration model beyond 2000 impossible. 
Following this spike, the response tended to stabilize in the period from mid 2002 to early 
2003. The response changed slightly in mid 2003 (corresponding to the changes in 
response affecting the other lamps), and resumed a fairly stable response pattern until 
early to mid 2004, when two response spikes appear to have occurred (the poor sampling 
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of data during this period makes a definite statement more difficult). Later in 2004, a 
decrease in response on the order of 6% to 8% was observed (depending on band), with 
an indication that the response was returning to a previous level. The character of the 
post-bumper mode response can be seen more clearly in the expanded plot shown in 
Figure 6. 
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Figure 6. Band 1 [010] Bumper Mode IC Response. There is an apparent discontinuity in response 
beginning in mid 2002, after the transition to bumper mode. The response tended to remain fairly 
stable until mid 2003, when the same phenomenon that caused significant spikes affecting LS [001] 
and LS [100] caused a slight decrease. The response stabilized after this ‘upset’ event until early to 
mid 2004, when two response spikes may have occurred. These spikes were followed by a short 
period of relative stability and a significant drop in response observed later in 2004. 
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Lifetime IC Response in LS [010], Bands 5 and 7 

Figures 7 and 8 show the Band 5 lifetime and post bumper-mode IC responses, 
respectively. Interestingly, the response maximum in late 2003 is approximately 1 DN 
greater than the level observed for the early 2003 maximum. There does not appear to be 
any event in the Band 5 responses in the other lamp states, or in the Band 1-4 responses 
in this lamp state, that might account for this change; it is possible that this may be an 
artifact introduced by the outgassing correction. Due to the relatively sparse sampling, it 
is difficult to estimate the decrease in response for the scenes acquired later in 2004; for 
Band 7, it is on the order of 3%. 
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Figure 7. Band 5 [010] Lifetime IC Response. 



Figure 8. Band 5 [010) Bumper Mode IC Response. The outgassing maxima for 2002 and early 2003 
are at a comparable DN level; the maximum for the cycle later in 2003 is approximately one count 
higher. Due to the poor sampling throughout 2004, it is somewhat difficult to estimate the overall 
change in response for the scenes acquired later in 2004; it would appear to be less than the Band 7 
estimate of approximately 3%. 


Conclusion 

Since the later-lifetime (post- 1998) changes in IC response do not appear to have 
any correlation to image or vicarious calibration data, it seems reasonable to conclude 
that they are due to changes occurring only within the IC subsystem. These changes 
suggest the beginning stages of failure either in the lamps themselves, or in the fiber optic 
cables and/or lenses that direct the lamp outputs to the focal planes. A cursory 
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examination of post-2002 telemetry data [1] suggests that at least through 2003, there did 
not appear to be significant changes in the lamp currents that would demonstrate lamp 
failure. It is possible that the lamp filaments are experiencing accelerated degradation, 
leading to random spikes/drops in their overall output as bits of the filament bum off. As 
LS [100]’s output is not attenuated, it would seem reasonable that any changes in its 
output would be easiest to detect; the attenuated nature of LS [001] and [010] could be 
‘masking’ significant changes to some extent. A more thorough sampling of IC data 
throughout the lifetime, as would be provided by the upcoming TM Image Assessment 
System (TMIAS) could help to answer some of these questions. 


References 

1. http://landsat-cal.gsfc.nasa.gov/ias web/L5 trending.html . A valid user ‘ID’ 
and password is required to access the site. I used ‘truggles’ and ‘tr_LPSO’ for the 
username and password, respectively. 
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Vicarious Calibration 

2002-2004 


The radiometric calibration phase of SDSU’s Landsat activities was a 
continuation of the work begun by Steve Schiller from the Physics Department. 

Beginning in 2002, this work was transitioned to a new team led by David Aaron (also in 
the SDSU PhysicsDepartment). The primary goal of the 2002-4 phase was to maintain 
the ground based vicarious calibration assessment of Landsat 5 and 7. The secondary goal 
was to fully transition the vicarious calibration activities to the new team. The final goal 
was to ‘productionize’ data collection and analysis procedures and processes. 

The vicarious calibrations were successfully accomplished; however, the recent 
cool cloudy summers in South Dakota minimized the number of clean ‘high confidence’ 
data collections. The project has been successfully transitioned to the new team at SDSU. 
This transition was accomplished with the kind assistance of Steve Schiller, the Goddard 
team, Kurt Thome’s University of Arizona Remote Sensing Group, and the JACIE group 
from Stennis Space Center. Finally, excellent strides were made in reworking our 
procedures, instrumentation set, and especially in upgrading the analysis software. 
Protocols now have the team routinely doing data verification in less than 24 hours of a 
collection and, while not yet routine, pilot studies have verified a process capability of 
less than 1 week for full data analysis. 

Activities: For 2002-2004 data collections were attempted for most overpasses during 
the summer months (June thru mid September) for both Landsat 5 and 7. Data collection 
activities centered around Brookings South Dakota. In 2003 an attempt was made to 
expand the data collection to a crop site (com). However, a pullback to the traditional 
‘scene center’ grass site was made as a response to the Landsat 7 SLC malfunction. 

Personnel from the radiometry team were included on trips to Goddard Space 
Center and NASA headquarters for technical interchange meetings. Also David Aaron 
and Gyanesh Chander (EROS Data Center), traveled to the University of Arizona for 
instrument calibration/repair and participated in a Ivanpah Playa data collection. 
Numerous meetings and joint data collections (primarily for high resolution commercial 
satellite assessment) with Stennis Space Center personnel were also undertaken. 

Radiometry; 

Data Sites and Equipment: 

Grass Site: 

The primary ground site used for radiometry in this study is within an area located in 
southeast Brookings that has been previously designated as ‘the 3M site’. The site is 
located within a field consisting of rough grasses and weeds. The site and the surrounding 
40+ acres contain the same vegetative mix to minimize adjacency effects. In previous 
years, radiometric irregularities due to uneven live vegetation and clumped dead 
vegetation have seriously complicated data analysis. To minimize this problem, an area 
within the site approximately 250m EW by 150m NS in 2002 and slightly reduced 2003- 
4 was predefined. To simplify MTF measurements (for related high resolution satellite 
work) the area is rotated from true north by 6 degrees. Maximum elevation change is 4.9 
m (see contour map). The site location (courtesy Stennis GRIT): 


12 



Northwest comer at: Lat: 44°17'31.12"N Long: 96°45'59.34"W 
Southeast corner at: Lat: 44°17'25.15"N Long: 96°45’48.70"W 


This area was routinely maintained by either a mowing service (2002) or by members of 
the calibration group. The ‘grass’ was nominally kept mowed to a height of about 10cm 
using tractor-powered 6-foot rotary mowers. In 2002 and for the August thru October 
2004 dates, the east 80m and the west 20m of the site was mowed using a ‘lawn mower’ 
and the clippings were bagged and removed. These areas were used for pulse MTF (‘blue 
tarp’) and mirror MTF studies for ancillary high resolution satellite projects. Within the 
center section the clippings were not collected. For safety, within the tarp areas an 
attempt was made to fill various in-ground holes (gopher/badger dens). 

To better maintain data reproducibility and facilitate setup, a system of in-ground 
markers (30 cm iron spikes labeled with plastic ‘ear tags’) was established and surveyed 
(courtesy of the Stennis GRIT). The markers were deployed to define the site boundaries 
and the location of select equipment and targets. These markers were driven flush with 
the ground (recall the site must be mowed). If a stake location was ‘lost’, the stake was 
easily re-located using a standard metal detector. A basic contour map of the site was 
developed in 2003. The availability of the site, owned by the Brookings division of the 
3M Corporation, is gratefully acknowledged. 

Instrumentation: 

The basic electronic/optical instrumentation used by the SDSU group has been 
previously described. The primary equipment base includes: 

• Automated SunPhotometers designed and built by the University of Arizona’s 
Electrical and Computer Engineering Department. 

o Unit 08: A good, solid instrument which has been in use for many years. It 
generally works well but is beginning to show its age (occasional heater 
malfunction and tends to lose tracking in moderate cloud cover). Always 
deployed on CEH rooftop (about 2 miles from site). 

o Unit 30: A new ‘2000’ series sunphotometer received in April of 2004. 
More compact and tracks well even during moderate cloud cover. The unit 
is typically field deployed using a generator as a power source. 

o Unit 20: This was borrowed from UA for part of the 2003 season (when 
unit 08 was experiencing problems) 

• A Fieldspec model FR spectro-radiometer (ASD unit 638). Work in 2002 was 
under calibration 8 and 2003-4 was under calibration 10. 

o Backup was generally provided courtesy of EROS Data Center and their 
unit 6219. The backup was never ‘required’; however, on some collections 
cross validation was done using unit 6219. 

• An 18 inch 99% reflectance Spectralon" 11 panel (Labsphere) is used for ‘white 
reference’ data. BRDF was characterized by the University of Arizona’s remote 
sensing group. 

o For backup, new 12” 99% and 12” 75% reflectance panels were purchased 
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• Several MFR Shadowband units (by Yankee Environmental Systems). Generally 
one is deployed at the rooftop site and another is deployed at the grass site 

• RTRM (Real time reference monitor). A new instrument was designed and built 
based around an Ocean Pacific compact spectrometer. This unit provides a 
relative but continuous monitor of radiance as reflected from a 75% spectralon 
panel. This data is used to then sort and interpolate the actual ASD ground data. 

• Ancillary support instrumentation for weather, sky path radiance, etc. 


In addition to the above, a computer with a backup RAID unit was procured and an 
archiving plan for field data was established. Data is also routinely backed-up to an on- 
campus mainframe archival system. 

Data collection Protocols: 

For the ASD grass data, the 8° optic was mounted on an adjustable pole and held 
at a nominal height 1.5 m above the surface (at normal to surface). This produces a 
sampling area of about 20 cm in diameter. For the white reference data, the panel is 
leveled and the seam is aligned with the sun to eliminate shadowing. The ASD optic is 
held about 50cm above the panel. Grass data is taken while walking on NS transects of 
the site. Generally from 8 to 10 ‘rows’ are collected, interspersing white reference panel 
data between each pair of rows. 

Procedure: 

-Warm up instrument for 15 minutes. 

- Collect ‘white reference’ data (generally three sets) 

-While walking, collect grass data down (south) one row and up (north) the next. 

ASD setting: 

Set to 20 spectra/file; 10 files per save 

(Typical Folders: Aug23\SDSU3M) 

This parameter set results in 40 or 50 files per row (depending on walk speed), 
-collect another white reference data series 
-walk/collect on the next grass rows, down and back 
-continue 

A complete collection over the site requires about 14 hour and produces 500-600 files. If 
there are ‘scattered’ clouds, collection is generally interrupted any time a cloud obstructs 
the sun. During data collection, a ‘walk sheet’ is generated (by hand), which attempts to 
describe what was happening and when (note that there may be minor discrepancies in 
file numbers on this sheet. Usually they are quite obvious). 

The 10 band sunphotometer (ASR) data is collected every minute from sunrise to sunset. 
Data collected with unit 08 is from the rooftop of a building (elevation 518m above sea 
level versus the site’s elevation of 505m). Starting in 2004, the new ASR 2000 (unit 30) 
is generally deployed at the north edge of the site. Shadowband MFR (for evaluating 
diffuse to global radiance ratio) units are always deployed at rooftop and sometimes at 
the north edge of the field. Vo values (zero air mass DN) for the sunphotometers are 
obtained and monitored by deploying the units on ‘clean’ days. For all units, stability has 
been excellent (generally <1% change per year for ASR 30, over 3 years for ASR 08). 
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Weather data is primarily from the Brookings municipal airport (~2km from site, 
elevation 502m). In addition, weather satellite imagery (GOES) is downloaded and 
archived for reference purposes. 

Two to three, blue ‘comer marker’ tarps are deployed just outside of the NW and the SE 
boundary to assist in clearly delineating the site area on the imagery. 

Data Processing: 

The raw grass spectral data is processed using a data reduction program which produces 
time and white reference panel BRDF corrected reflectances averaged over the entire site. 
This is a MATLAB based program that was developed and refined by SDSU personnel 
(2001-2003). The raw data is stored in the archive as ASD Level O and the reflectance 
data is stored as the ASD Level l. 

Sunphotometer data is reduced via the Langley analysis method. Again a MATLAB 
based program is used (developed 2002-2004). A front end is nearing completion that 
will semi-automatically ‘remove’ clouds (simply eliminating data lines where the sun 
was obscured by scattered clouds). 

The reduced sunphotometer data is passed to another program (presently an Excel 
spreadsheet, but will be folded into the MATLAB Langley analysis) which removes 
Rayleigh scattering and produces extinction coefficients for the 10 sunphotometer 
wavelengths. 

With the calibration activity being transitioned to new personnel, the group’s experience 
in MODTRAN modeling was severely degraded. Consequently, the group would like to 
thank and commend Kurt Thome and the Remote Sensing Group of the University of 
Arizona for their work in doing the modeling for the 2001-2003 data collections. The 
project phase simply could not have been continued without their assistance. 

Starting in late 2003 with full implementation in 2004, the SDSU group once again has 
atmospheric analysis capability. Primary analysis is under MODTRAN 4, with 
supplementary analysis (within an experiment to produce a very short turnaround, 
moderate accuracy model) using an Excel spreadsheet based model from NREL. 
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Figure 1 : A typical view of the SDSU grass site area from Sept 7, 2002 (Quickbird imagery). The area 
within the dotted lines is approximately 250 m east-west by 150 m north-south and is regularly mowed. 
Note the discrete targets shown are not generally deployed for Landsat acquisitions. 



Figure 2: Site elevations in meters as measured by SSC GRIT. Scale is approximately equivalent to that 
shown in figure 1. 
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Figure 3: Average grass hyperspectral reflectances for select dates in 2002. Note that the major water 
bands have been removed for clarity. 


Summary of SDSU Landsat Data Collects: 

Summer 2002 

Generally 12:04 & 12:05pmCDT overpass 
June 20: Landsat 7 & EO-1 
Good collect 
Light Cirrus Wisps 

Deployed Blue MTF tarps (Ikonos width) 

Ground ASD Reflectance Uniformity of: 6.2% 450-1350nm 

6.9% 1410-1800nm 
10% 1940-2470nm 

July 22: Landsat 7; EO-1 and Ikonos 
Good collect 

Good weather, slightly hazy 
Stennis Tarps Deployed 
Blue MTF deployed 
Mirrors Deployed 

Also deployed (N of maintained area), Plant Science ‘4 reflectance tarps’ 
ASD data from Stennis also 
Also extensive Cropscan II data 

Sept 8: Landsat 7; EO-1 (note a good Quickbird collect on the previous day) 
Good collect 

Reasonable weather, but hot & humid so somewhat hazy 
Only corner marker tarps deployed 

Recorded ancillary data on spots usually covered by Stennis & MTF tarps 
Ground ASD Reflectance Uniformity of: 5.4 % 450-1350nm 

6.3% 1410-1800nm 
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14.3% 1940-2470nm 


Summer 2003 

Landsat 5: 2003 Useable Data Collections 

May 30: ‘Brookings Site’ Good but very windy (dust at ground level) 

Sept 3: 3M Good 

Sept 19: 3M Very good 


Landsat 7: 2003 Useable Data Collections 

July 25: 3M Fair (hot and hazy) 

Aug 26: 3M Good 


Summer 2004 

Landsat 7: 2004 Useable Data Collections 

June 25: 3M ‘shoot between clouds’ comparison with EOl 

Fair+ (“good” for relative quantification) 

S. Schiller did initial analysis using ‘light trapping’ 
model 

July 1 1 : 3M Fair : scattered clouds/ cirrus over site 

Sept 29: 3M Very good 


Gain Calculation Summary: 
Landsat 5 



2003 

3-Sep 

19-Sep 

Gain 

*variation 

Gain 

* variation 


0°) 


(la) 

L5-B1 

1.295 

2.1% 

1.285 

3.6% 

L5-B2 

0.653 

3.2% 

0.667 

3.8% 

L5-B3 

0.901 

2.2% 

0.886 

4.3% 

L5-B4 

1.051 

1.9% 

1.046 

2.1% 

L5-B5 

7.994 

6.7% 

7.293 

6.6% 

L5-B7 

13.681 

5.6% 

13.362 

7.3% 


L5-Gain is in DN per W/m2 sr pm 
*NOTE: Listed variation is image 
based only 
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Landsat 7: ETM+ Gains 2002-2004 SDSU 
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Conclusions: 

The trending for both the TM and the ETM+ sensor systems are consistent with 
previous calibration activities. No change in the present gain coefficients is 
warranted. 

Data acquisition and analysis techniques at South Dakota State University have been 
fully restored and procedures are now streamlined to improve turnaround time and to 
make the process more robust. 

In conjunction with continued data collection and analysis for the grass site, two 
additional tasks have been initiated. The first task will be to develop a formal 
assessment of experimental uncertainty for the SDSU data collection/analysis. The 
second task is the refinement of a formal algorithm for data analysis in the presence 
of scattered cloud cover. This will include implementation of the data sorting 
algorithm mentioned previously and in continuing development of a ‘light trapping’ 
model to accommodate the changes in diffuse radiance due to scattered clouds 
reflecting upwelling radiance. 
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Relative Gain Characterization 


1.1 Background 

Relative Radiometric Gain (RG) refers to the phenomenon that occurs whenever 
an imaging system has more than one detector. With multiple detectors, the responsivity 
of each detector is slightly different from its neighbors, and possibly every other detector 
in the array. A model that is often used is 

DN = gj*r + bj (1) 

where DN refers to the digital number recorded in the image, gj is the gain of the ith 
detector (in DN/radiance), r is the radiance incident on the detector, and bj is the bias of 
the ith detector (expressed in units of DN). Thus, if a number of detectors are exposed to 
a uniform radiance level, r, the response will not be uniform. As a result, significant 
striping visually exists in the imagery and is most pronounced in the homogeneous 
regions. Since the bias is recorded after each scan, it is a known quantity for each detector 
that can be removed from the equation and only gj needs to be estimated for each detector 
[ 11 - 

In this document, a new methodology to destripe Landsat 5 scenes is discussed. 
The objective is to develop a definitive characterization of the relative radiometric gain of 
all reflective channels of the Landsat 5 Thematic Mapper instrument. The goal was to 
develop quantitative models for each of the detectors/channels in the Thematic Mapper. 
A scene-based approach is used in order to calculate relative gain for each detector. 
Relative gain (RG) is calculated as (Detector mean / Scene mean). The scenes used in this 
analysis were artifact corrected before the relative gain calculations were performed. The 
relative gain values are then modeled into Oth order or 1st order linear models for all 
detectors/channels. These statistical models are then applied to all the detectors and are 
used to correct the scene data. The relative gain correction removes residual striping 
present in scene data. 

1.2 Relative Gain Calculation 

The underlying theory behind histogram equalization is that the image data from 
all detectors should have the same mean and standard deviation if they are exposed to the 
same radiance source. This realization gives the foundation for balancing the radiometric 
response of each detector [2]. The following development forms the basis for RG 
correction by histogram equalization. Assume DNj and DN2 are the digital output values 
of two detectors, detector 2 being the reference detector. 

Then 

DNi = gir + bi and DN 2 = g 2 r + b 2 (2) 

If both detectors observe the same radiance, then the output DN values should be equal. 
Substituting for radiance from the first equation into the second yields 

DN 2 = (g 2 /g,)x(DNi-bi)+b 2 (3) 

which is the relationship used to normalize the output of detector 1 to the reference 
detector. The mean value from each detector is 

DNimean = r mean gi + bi and DN 2mean = r me ang2 + b 2 (4) 
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It can also be shown that the ratio of the gains is equal to the ratio of the standard 
deviations [2], 

ai la 2 = g t /g 2 (5) 

Thus, to implement relative gain correction via histogram equalization, the ratio 
of standard deviations is substituted into (3) and, in the case of Landsat TM, bias is 
estimated from the calibration file associated with each scene. However, since bias is 
accurately estimated from the calibration file, another approach exists to determine the 
ratio of relative gains. Consider first removing the detector bias to obtain a net DN value 
as follows: 

DNi’ = DNi -b, = g,r 

and DN 2 ’ = DN 2 - b 2 = g 2 r (6) 

Then the means of the net DN value are simply 

DN lmean — glTmean and DN 2 mean — g 2 rmean 

which leads to 

gi/g 2 = DN’ lmean/DN 2mean ■ 

Thus, the gain ratio can also be calculated from the ratio of means of the net DN values. 
The lifetime relative gain characterization was implemented using this method where the 
reference detector was the mean of all detectors (i.e. the scene mean). This approach was 
felt to be less sensitive to noise since first moments are used rather than second moments. 

1.3 RG Characterization 

Before computing the relative gain at a scene level for individual detectors, the 
unprocessed Level 0 (LOR) data were subject to Scan-Correlated Shift (SCS) and 
Memory Effect (ME) artifact corrections [3] using the processing system implemented at 
South Dakota State University. Detector bias on a line-by-line basis was calculated from 
the artifact-corrected calibration data through averaging over a 550-pixel window in the 
shutter region; these bias values were then subtracted from the artifact-corrected image 
data. Because Memory Effect had been corrected previous to this step, the bias estimate 
should be accurate for the entire scan [4], 

The calculation of detector and overall scene means used in the relative gain 
estimates relied on a restriction of image data to values between 5 DN and 245 DN. This 
was done in an attempt to remove from consideration ‘saturated’ data and/or potential 
negative radiance values due to bias error. Although no explicit cloud mask was used, the 
restricted DN range tended to remove clouds and shadows, which helped in stabilizing 
the analysis. 

From the detector and scene statistics for each scene, the relative gains for each 
detector were calculated. These time series were used to develop linear gain models for 
each detector as follows: 

• The relative values were first plotted versus time in units of days since launch 

(DSL), as shown in the example of Figure 1. 

• The data points were then subjected to outlier detection. Those with values 

outside the range of g m eani(t) ± 2*a were considered as outliers (where gmeani(t) is 
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Relative Gain Relative Gain 


the relative gain value of the ith detector and a is the standard deviation of the 
corresponding data values). 

• First order linear functions of the form gj(t) = An*t + A2i (where Ah and A* are 
the slope and intercepts of the ith detector) were fit to the data. 


YUMA, Band 3, Detectors 1-4 YUMA, Band 3, Detectors 5-8 



Days Si nee Launch Days Since Launch 

(a) (b) 


YUMA. Band 3. Detectors 9-12 YUMA ' Band 3 > Detectors 13- 16 



Days Since Launch Days Since Launch 

(c) (d) 


Fig I Linear Model for Detectors in Band 3 (Primary Focal Plane), (a) Detectors 1-4 (b) Detectors 5-8. (c) Detectors 9-12. (d) Detectors 13-16. 
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1.4 Validation of Relative Gain Models 

Three scenes of the White Sands, NM test site (WRS path 33, row 37) acquired 
May 11, 1986, April 20, 1990, and June 1, 1998 were chosen for correction using the 
proposed relative gain models. Each validation scene was processed in the same manner 
as the scenes used in the modeling procedure. The scenes were then corrected using date- 
specific relative gains predicted from the models. For purposes of comparison, the 
scenes were also corrected using standard histogram equalization. Figures 2a, 3a, 4a 
show a 400x400 pixel area extracted from the common subregion for each scene. In 
these images, striping is clearly evident. Figures 2b, 3b, 4b show the same area after 
correction with relative gains estimated from the lifetime models, and Figure 2c, 3c, 4c 
show the same area corrected with histogram equalization. All images have been 
identically stretched before and after correction to emphasize potential changes resulting 
from the correction. In general, the images corrected with relative gains derived from the 
lifetime models show less evidence of residual striping than do the images corrected with 
histogram equalization. 


To provide a quantitative estimate of the degree of residual striping present after 
correction, an ‘integrated’ striping metric was calculated for a 400 x 400 pixel 
homogeneous area within the common subregions before and after correction: 


SM = 


ISR-before ~ after 


ISR, 


1 before 


x 100% 


( 7 ) 


The ISR terms in (7) represent an ‘integrated’ striping ratio, defined as the mean striping 
ratio of all ‘significant’ components found at multiples of the fundamental frequency of 
1/16 cycles per pixel (up to and including the Nyquist frequency): 



(8a) 


The striping ratio components are estimated from the image subregion, and are defined as 
the ratio of the average Fourier transform magnitude in the along-track (across-scan) 
direction to the average Fourier transform magnitude in the cross-track (along-scan) 
direction [5]: 


SR n 


|f(0| 


(8b) 


An example striping ratio spectrum before and after relative gain correction is shown in 
Figure 5. Tables I through VI show the estimated ISR before correction, after relative 
gain correction, and after histogram equalization; lower values indicate a greater degree 
of correction. Overall, differences between the two methods tend to be fairly small. 
Performance of the new approach was better in the primary focal plane bands than in the 
cold focal plane bands, consistent with the adequacy of the modeling. However, 
implementation of the relative gain approach simply involves use of a lookup table, 
whereas histogram equalization requires calculation of image statistics. 
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Stripe Ratio 



(a) (b) (c) 

Fig 4. WhiteSands Scene (June 1998) - Band 3 (a) Before (b) After RG Adjust, (c) After Histogram Equalization 



Spatial frequency- cycles/sample Spatial frequency- cycles/sample 

Fig.5. Striping Ratio of WhiteSands Scene (June 1998) - Band 3 (a) Before RG Adjust, (b) After RG Adjust. 
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TABLE I 

ISR ESTIMATES FOR BAND 1 


Scene 

ISR Before 

ISR After RG 

ISR After Hist. 

SM 


Location/Date 

Adjust 

Eq. 

RG adjust 

Hist. Eq. 

WS- 86,131 

6.001 

0.013 

0.272 

99.778 

95.469 

WS-90,110 

5.436 

0.277 

0.098 

94.901 

98.197 

WS- 98,180 

8.839 

1.691 

2.375 

80.865 

73.131 



TABLE II 






ISR ESTIMATES FOR BAND 2 



Scene 

ISR Before 

ISR After RG 

ISR After Hist. 

SM 


Location/Date 

Adjust 

Eq. 

RG adjust 

Hist. Eq. 

WS - 86,131 

4.769 

0.341 

0.382 

92.852 

91.987 


2.352 

0.175 

0.221 

92.574 

90.626 


8.262 

2.049 

1.782 

75.196 

78.432 



TABLE III 






ISR ESTIMATES FOR BAND 3 



Scene 

ISR Before 

ISR After RG 

ISR After Hist. 

SM 


Location/Date 

Adjust 

Eq. 

RG adjust 

Hist. Eq. 

WS- 86,131 

12.990 

0.899 

1.030 

93.078 

92.069 

WS-90,110 

13.492 

1.458 

1.340 

89.191 

90.063 

WS- 98,180 

21.161 

0.639 

2.273 

96.983 

89.256 



TABLE IV 






ISR ESTIMATES FOR BAND 4 



Scene 

ISR Before 

ISR After RG 

ISR After Hist. 

SM 


Location/Date 

Adjust 

Eq. 

RG adjust 

Hist. Eq. 

WS- 86,131 

4.4475 

0.6664 

0.4922 

85.016 

88.933 


4.8312 

1.0436 

0.7458 

78.399 

84.562 


9.6472 

2.9768 

1.9643 

69.143 

79.638 



TABLE V 






ISR ESTIMATES FOR BAND 5 



Scene 

ISR Before 

ISR After RG 

ISR After Hist. 

SM 


Location/Date 

Adjust 

Eq. 

RG adjust 

Hist. Eq. 

WS- 86,131 

7.223 

0.551 


92.369 

97.701 

WS-90,110 

7.001 

0.976 

0.424 

86.064 

93.950 

WS- 98,180 

21.748 

2.488 

1.177 

88.559 

94.590 



TABLE VI 






ISR ESTIMATES FOR BAND 7 



Scene 

ISR Before 

ISR After RG 

ISR After Hist. 

SM 


Location/Date 

Adjust 

Eq. 

RG adjust 

Hist. Eq. 

WS- 86,131 

7.115 

0.575 


91.921 

99.525 

WS-90,110 

6.202 

0.802 

0.131 

87.072 

97.882 

WS - 98,180 

6.658 

2.138 

1.885 

67.884 

71.684 


28 










1.5 Conclusions 

The analyses presented in this report demonstrate the stability of the Landsat 5 
Thematic Mapper over a 20-year period of continuous operation. Relative gain 
differences between detectors vary from approximately 0.97 to 1.02. These gains, with a 
few exceptions, have either been constant or slowly varying over the lifetime. Those 
detectors that have changed slightly have done so in a manner that allows 
characterization through fitting of simple linear models. However, detector 1 in Band 7 
could not be properly characterized through a linear model. These results have led to an 
improved method of relative gain correction that is based only on the date that the image 
is acquired instead of scene content. The RG method is also computationally much 
simpler than traditional histogram equalization. 

However, there are detectors in bands 5 and 7 that have a noisier response thereby 
decreasing the modeling accuracy. Hence, this would impact the degree of relative gain 
correction that could be achieved when applied to the scene data. Quantitative results of 
Band 4 suggest that histogram equalization always tended to produce a higher degree of 
correction over the suggested lifetime approach. Based on the quantitative results shown, 
it can be concluded that the lifetime model approach produced a higher degree of 
correction approaching 75% to 99.7% in Bands 1, 2, and 3. For Bands 4, 5, and 7, 
application of the same approach results in removal of 67% to 92% of striping. With this 
approach, relative gain correction can be more accurately performed and overall 
radiometric fidelity of image data acquired by Landsat 5 Thematic Mapper can be 
preserved. 
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Outgassing 


Problem description 

An artifact, commonly known as Outgassing plagues the radiometric responsivity 
of the TM Cold Focal Plane (CFP) bands. Detector responses to Internal Calibrator (IC) 
pulses in Landsat 5 are observed to fluctuate in a cyclic manner about 4.5% in Band 5, 
and about 3.5 % in Band 7, and exponentially decay in Band 6 over single outgassing 
cycles - periods between two outgassing events - throughout the instrument's lifetime. 
Figure 1 represents an example of Bands 5 and 7 detector responses observed during an 
outgassing cycle in 1987. The attenuation in detector responsivity is removed and every 
detector’s performance regained whenever an outgassing procedure - warming the 
housing of the CFP detectors from its normal operating temperature of -1 81° C to the 
temperature of the primary focal plane (approximately 23°C) followed by cooling down 
with the purpose to remove potential contaminants within the CFP optics - is performed. 
Thus, it is postulated that a thin film of some material is deposited on the window 
surfaces of the CFP dewar much like frost builds up on automobile windshields during 
winter weather. 


Outgaaaing Data 1987. Band 5 


1120 1140 1180 1180 1200 1220 1240 

Tima, DSL 


Outgaaaing Data 1987, Band 7 


42.8 o o 


Vi 20 1140 1180 1180 1200 1220 1240 

Tima, DSL 


Fig. 1. Bands 5 and 7 detector responses observed during an 

outgassing cycle in 1987 


Theory 

Two-film interference model is developed to simulate the cyclic pattern present in 
the mid-infrared detector responses of Landsat 5 TM. This model is able to predict Band 
5 and 7 responses at any point in time simply by knowing when the last outgassing event 
occurred. 


30 


The two-film model layout representing the contaminated ZnSe CFP dewar 
window is depicted in Figure 2. An antireflective coating (ARC) deposited on the ZnSe 
substrate to increase transmittance through the window is included in this model as the 
first thin film layer. The second layer represents the unknown contaminant. An 
electromagnetic wave incident on such a structure is partly reflected at each of the 
boundaries, I, II, and III. Additionally, the wave is reflected backward and forward 
between the interfaces. Multiple reflected waves combine with the incident waves in a 
process of interference. 


Vacuum 

yo, n 0 

Contamination 

yi, Ni, di,S, 

ARC 

y 2 , n 2 , d 

ZnSe 

i /\ /\ /\ /\ /" 

y3, n 3 
/\ /\ / 


I 

II 

III 


Fig. 2. Two-film model layout. The boundaries between 
layers are labeled with Roman numerals. The ZnSe layer 
represents the dewar window material. 


The mathematical treatment used to derive expressions for electric and magnetic 
fields at the interfaces is based on an application of Maxwell equations at the interfacial 
boundaries. For a normally incident electromagnetic wave at interface I, the electric ( E ) 
and magnetic ( H) fields on the corresponding boundaries I, II, and III are related as: 


E, 

H, 


COS 5 X 

sin£. 

J 

y\ 

Jy\ sin<J i 

cos 5 X 

- 

sin S , 

cos 5 2 

J 

y 2 

_jy 2 sin S 2 

cos S 2 


y 2 E n 


= M. 


“ II 

yi E n 



f 

^ III 

= m 2 

E m 


yi E m_ 


y 2 E ui. 


( 1 ) 

( 2 ) 


where Mi is the characteristic matrix of the film buildup that has the complex index of 
refraction N x -n x - jk x having real index of refraction «/ and imaginary part, ki, that 
accounts for absorption by film material (index of absorption), and the characteristic 
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admittance y, = N l Y 0 , where Yo = 1/377 Siemens is the admittance of free space. M 2 is 
the characteristic matrix of the ARC that has the characteristic admittance y 2 = n 2 Y 0 and 

is assumed to have a real index of refraction ri 2 , and ys is the characteristic admittance of 
the ZnSe substrate in which the absorption is neglected (ns is real). The phase difference 
due to wave propagation through the film build up is given by: 

s _ 2 gjVi ^ M”. - AM (3 

1 A 2 1 

and the phase difference that occurs when the wave traverses the ARC is: 


_ 2 m 2 d 2 
O') — ■ 


A 


(4) 


where d t and d 2 stand for the thickness of the contamination and ARC layers, 
respectively, and A is the operating wavelength. 

The relation among the fields on the interfaces III and / is then: 



= M, 

' E a ' 

= m.m 2 

’ E,„ 

= M 

F 

^ 111 



1 

h? 

rs 


_y 3 e ]jj _ 


1 

t*} 


(5) 


where M is the characteristic matrix of the two-film model: 


M = 


y 1 1 

cos d i cos S 2 — —sin <5, sin <5 2 j — sin S 2 cos <5,+ j — sin£,cos£ 2 


Ti 


y 2 


y, 


Jy, sin d, cos d 2 + jy 2 cos <5, sin£ 2 sin d, sin d 2 + cos 5 { cosS 2 

y 2 


( 6 ) 


After normalizing to the electrical field at the rear interface, Em, the fraction of 
the irradiance, T, transmitted to the ZnSe substrate can be derived from boundary 
conditions as: 


4 To Re(y 3 ) 

(y 0 E ln +H ln )(y 0 E In+ H ln )* 


(7) 


where and Hj n are the normalized electric and magnetic fields at the front interface, 
respectively and represents the complex-conjugate operation. Once the fraction of the 
transmitted power is determined, it has to be multiplied by a scaling factor specific for a 
given IC Lamp State (LS), Sf to predict the detector responses, DN: 

DN = Sf xT (8) 
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Outgassing Modeling 

Four independent, closely-sampled scene sets representing outgassing cycles in 
years 1985, 1987, 1993, and 1999 have been acquired and the average Net-Integrated 
Pulse Values (NPV) in Digital Numbers (DN) derived from odd numbered detector 
responses to the second strongest IC lamp (LS 010) during reverse scans. The two-film 
model parameters are estimated based on analysis of each cycle separately in order to 
achieve consistent model parameter values and, thus, to allow usage of the model for 
lifetime characterization. 

The contamination material, as well as the other media employed in modeling of 
interference effects on the CFP, was assumed to be linear, homogenous, isotropic, and 
nondispersive across the bands. Furthermore, the operating wavelengths are assumed to 
be equal to the corresponding band center wavelengths. That is: 1650 nm for Band 5 and 
221 5 nm for Band 7. The refractive index of vacuum (free space) is known and equal to 
1 for all wavelengths. The refractive index of Zinc-Selenide is also known and equal to 
2.45 for Band 5 and 2.44 for Band 7. The other model parameters have to be estimated 
using a nonlinear curve fitting technique. One of the keys to success in such a process, 
which is by nature an iterative one, is to define sufficiently good starting parameter 
values. 


The thickness of the film buildup has to be expressed as a function of time. From 
the nature of the interference phenomenon, it is known that the optical thickness of the 
film buildup, nidi, increases for Ml in a time period between two consecutive local 
minima of transmitted energy. This explains the different periods of oscillations in Band 
5 and 7 responses (Figure 1). This also means that there is a direct relationship between 
the film growth rate and the periodicity of the detector responses. In order to characterize 
changes in period over the lifetime, the initial periods were estimated by minimum- 
squared-error fitting of exponentially attenuated cosine curves to the outgassing data sets. 
The oscillation period was not found to change significantly between the 1985 and 1987 
data sets. Later in life (post- 1987), a reasonably linear increase of period (and a 
consequent decrease of the film buildup rate) was indicated in both bands, which can be 
expressed as a function of DSL: 


T p (/) = mt + b (9) 

where m and b are the slope and intercept, respectively, of a regression line fitted to the 
estimated oscillation periods, and t is the time since launch in DSL. Thus, the initial 
parameter values for m and b in Band 5 are 0.03968 and -4.823, and in Band 7 0.06167 
and -20.62. 


The film thickness will increase during a given period, T p , by: 



( 10 ) 
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The film buildup rate, n, in nm/day is the ratio between the given film thickness 
cf^and the time needed for its deposition and, is given by: 


r *(0 = - 


l 


( 11 ) 


for the early lifetime (T p constant), and 


r b (0 = 


'in. 


T p (t) mt + b 

for the later lifetime, taking into account the time dependence given by (9). 


( 12 ) 


The film thickness, d/(t), accumulated since an outgassing event occurred on day 
t 0 G can be determined by integrating the expression for the film build-up rate. For the 
early lifetime that is: 


d x (0 


2 n x T p 



and later in life: 


d x (t) 


1 

2n x m 


f 

In 

V 


mt + b 
rntoo +b 




J 


(13) 


(14) 


The contaminating film material is not known, but there is some evidence that it 
could be ice or an ice mixture that forms by condensing water vapor on the Cold Focal 
Plane window. The wavelength dependent complex refractive index of ice for the center 
wavelength of Band 5 (1650 nm) is 1.2878-jO. 000241, and for the center wavelength of 
Band 7 (2215 nm) it is about 1.2608-j0. 000232. 

Information about the structure of the ARC used for coating of the CFP ZnSe 
windows in Landsat TM was not available. For that reason, the initial values of the 
effective thickness and refractive index were estimated according to general ARC design 
principles — zero reflectance is achieved when the ARC has an index of refraction equal 
to the geometric mean of the refractive indices of the adjacent media, and an optical 
thickness (product of the thickness and the index of refraction) equal to an odd multiple 
of one quarter of the operating wavelength — and with regard to the band center 
wavelengths. Then, the initial values for Band 5 are n ARC = ^n 0 n /nSe - f\ ■ 2.45 =1.565 
and d 2 = 263.6 nm, and for Band 7 n 2 = 1.562 and d 2 = 354.5 nm. 

The value for the scaling factor, S </, used in Equation (3.8) needs to be close to the 
average detector responses in DN since the transmittance, T, is always close to 1 . Then, 
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the appropriate initial values for LS 010 would be in the range of 32 - 35.5 for Band 5 
and 42.5 - 47.5 for Band 7. 

The final parameter values for each single outgassing data set were estimated 
based on a least-squares error criterion that is optimized using the Nelder-Mead downhill 
simplex algorithm. There were many sets of parameters for each outgassing data set that 
provided very similar solutions using a least-squares error criterion. Therefore, as 
additional condition, the final parameter set is to be chosen such that the contaminant’s 
index of refraction is closest to the index for ice. 

In order to extend procedures used in modeling of single outgassing cycles to the 
whole instrument’s lifetime, the outgassing data sets had to be corrected (normalized) for 
the apparent increase in overall Band 5 and 7 detector responses after approximately 
1600 DSL due to changes in IC lamp outputs. The correction for this Tamp effect’ (LE) 
is illustrated in Figure 3. First, mean values for each data set were calculated from the 
responses acquired during the time between the outgassing procedure and the following 
maximum, and a quadratic function,^, was defined such that it passed through those 
mean values. The 1993 and 1999 data sets were then normalized to the 1987 data 
through multiplication by a correction factor, Clt, defined as: 


C lt (DSL) 


A*S7 

fu(DSL) 


(15) 


where Ms 7 is the mean DN value of the first cycle in 1987 response data. After the LE 
correction, the variability in Band 5 and 7 data is assumed to be only due to outgassing 
effects 


Lamp Effect Correction for Band 5 Odd Detectors 
35 
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Lamp Effect Correction for Band 7 Odd Detectors 
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TOO 2000 3000 4000 3000 6000 

Days Since Launch (DSL) 


Fig. 3. LE correction function fti applied to Band 5 and Band 7 NPV data 


Fitting the two-film model to detector response data from individual outgassing 
cycles resulted in slightly different parameter sets for each cycle. As the oscillatory 
pattern in the Band 5 and 7 responses is assumed to be caused by accumulation of the 
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same contaminant over the whole lifetime, model parameters for each outgassing cycle 
should remain constant. To find a unique parameter set applicable over the entire 
lifetime, a full lifetime model (LM) was constructed using all available individual data 
sets, with the oscillation period remaining constant early in life and allowed to increase 
later in life. DSL 1434 was chosen as the approximate date when the oscillation period 
began to increase, as it corresponds to the date imaging was resumed after outgassing 
procedure #18. The date when the detector responses start to increase due to changes in 
lamp output is also not precisely known. For the purposes of this analysis, DSL 1630 
was chosen as the last date prior to its onset. 


Modeling results 

Final parameter estimates derived for each individual outgassing data set and for 
the lifetime model fitted to all single cycle data sets together are shown in Tables 1 and II 
for Bands 5 and 7, respectively. The resulting Band 5 and Band 7 LM fits are shown in 
Fig. 4. In the lower four plots, the upper curve represents the original data, and the lower 
curve is the LE corrected data with the model fit. The overall mean squared error (MSE) 
of the lifetime model indicates that a better lifetime fit was obtained for Band 7. This is 
primarily due to the poorer fit in the early lifetime for Band 5. Not surprisingly, this is 
the time interval with the shortest oscillation period in general. In addition, the Band 5 
oscillation period is even shorter than the corresponding Band 7 period, thus it will be 
more sensitive to parameter estimation errors. 

Table I: Estimated Model Parameter Values Derived from 

Nonlinear Fits to the Single Outgassing Data Sets and LE- 
Corrected Lifetime Model (LM) , Band 5 


Paramete 

r 

1985 

1987 

1993 

(reduced) 

1999 

LM 

«/ 

1.2878 

1.2881 

1.2879 

1.2891 

1.2878 

k, 

0.000879 

3 

0.000166 

0 

0.00006549 

0.0010 

87 

0.00072 

58 

r >2 

1.683 

1.6420 

1.6759 

1.6849 

1.6739 

d 2 , nm 

272.1 

257.5 

272.2 

269.4 

269.0 

m 

- 

- 

0.04091 

0.0418 

8 

0.03876 

b 

- 

- 

-5.73 

-14.93 

-0.94 

T p , days 

43.70 

44.16 

- 

- 

45.75 

Sf 

32.382 

32.326 

34.018 

35.065 

32.41 

MSE 

0.003376 

0.005704 

0.006552 

0.0018 

38 

0.01787 


36 


Table II: Estimated Model Parameter Values Derived from 

Nonlinear Fits to the Single Outgassing Data Sets and LE- 
Corrected Lifetime Model (LM) , Band 7. 


Parameter 

1985 

1987 

1993 

(reduced) 

1999 

LM 

ni 

1.2608 

1.2606 

1.2603 

1.2606 

1.2606 

ki 

0.001768 

0.001646 

0.001530 

0.002493 0.002472 

n 2 

1.6733 

1.6789 

1.6658 

1.6633 

1.6677 

cl 2 , nm 

326.5 

311.9 

326.3 

331.5 

326.9 

m 

- 

- 

0.06221 

0.06207 

0.06224 

b 

- 

- 

- 20.93 

- 7.319 

- 18.59 

T p , days 

64.09 

63.44 

- 

- 

68.54 

Sf 

42.923 

42.977 

45.1036 

46.8801 

43.015 

MSE 

0.002616 

0.001795 

0.01386 

0.003350 

0.01111 


Lifetime Model Fit to Bend 6 1985 Outgassing Data Set 

32.4 



Lifetime Model Fit to Band 7 1986 Outgassing Data Set 
43 



Time, DSL 
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Ufttim* Model Fit to Bond 5 1993 Outgoing Data Sot 



3 °3050 3100 3150 3200 3250 3300 3350 3400 

Time. 0SL 




Fig. 4. Lifetime Model fits to individual outgassing data 
sets. The 1993 and 1999 data sets were corrected for LE 
prior to modeling (bottom curves and data points). In all 
cases, the MSE values are associated with the models. 


It is evident from the Tables I and II that the estimates for the effective optical 
properties of the ARC obtained from both models agree among all the data sets, with the 
exception of the 1987 data set; the difference is greater in Band 5. The ARC optical 
thickness can then be used to account for phase shifts present in detector response 
(provided that complete contaminant removal occurs during an outgassing procedure, and 
that the dates of outgassing procedures are accurately known). The different phase 
observed in the 1987 data set could be explained as a change in the ARC composition 
(e.g. deterioration of one or more layers in a multilayered structure), or by incomplete 
removal of the ice film during the outgassing procedure; however, the 1987 data set 
appears to represent the most complete outgassing procedure, as the detector response 
data start closest to a maximum. Differences in the absorption factor ki of the film 
thickness can be considered as a random error in modeling, because the estimation of its 
relatively small value is limited by the accuracy of the observed IPVs. 
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The slope, m, which characterizes change in oscillation period is found to be a 
very unstable parameter in each single cycle data set. It may change its sign with a very 
small change in the starting value of another parameter; this is usually followed by an 
abrupt change of the intercept b. This indicates that the change in period of oscillation 
that occurs during an OG cycle is so small that it is difficult to estimate well. Assuming 
that the amount of available film material decreases in an OG event (e.g. through a 
valve), there is a possibility that the period does not change at all during OG cycles. In 
that case, the observed change in period over the instrument lifetime could be explained 
as a step-wise change occurring after OG processes. 


Correction for Outgassing effects 

The DN value extracted from the lifetime model for the day of any outgassing 
event gives the detector response DNo, when no ice is accumulated on the CFP (32.1 DN 
for Band 5 and 42,83 DN for Band 7, when odd detector responses to IC lamp [010] 
analyzed). The DN value corrected for outgassing effects, DN cor , for any detector 
response ZW can be determined as: 

DN cor = DNx = DN x T ° = DN x— (16) 

c DN 0l0 SfxT T 

where DNoio is Net Integrated Pulse Value derived from odd numbered detector 
responses to IC lamp [010] at a given DSL, T 0 is the fraction of the EM radiation 
transmitted through the CFP without an ice layer, and T is the transmittance extracted 
from the LM for a given DSL. 


The 91 -scene SCS-corrected lifetime data set was used to test the corrections 
produced by the lifetime model. Results of the correction applied to the data set for 
Bands 5 and 7 are shown in Figure 5. From this figure, relatively good correction can be 
observed in both bands, with the exception of the early life period and the time interval 
during 1993 and 1994. Overall, the analysis results indicated that the variation from the 
estimated DNo value after correction was reduced 68% in Band 5 and 77% in Band 7 over 
the instrument lifetime after DSL 175, excluding the time interval during 1993 and 1994 
from consideration. Band 7 responses tended to be better characterized due to a longer 
oscillation period and consequent decreased sensitivity to phase shift effects. 
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Fig. 5. Lifetime Model correction applied to the Band 5 and 
Band 7 lifetime test data sets. 


Summary 

The correction for outgassing effects has been implemented for the Band 5 and 7 
absolute calibration lookup tables used by the National Landsat Archive Processing 
System (NLAPS) at the EROS Data Center. General application of this correction can 
potentially improve Band 5 and 7 calibration accuracy to levels consistent with the 
accuracy of the primary focal plane bands, on the order of 5%. 
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Landsat 4 Absolute Calibration 


A complete characterization and calibration was attempted for Landsat 4’s 
Thematic Mapper. This entailed analyzing the Internal Calibrator (IC) to verify the 
system was stable throughout the life of the instrument, applying existing correction 
models to the data to remove Scan Correlated Shift (SCS) and Memory Effect (ME) 
artifacts, and determining and applying a Relative Gain (RelGain) model to the data. 

Then, when the Landsat 4 data was completely corrected, a cross calibration model using 
Landsat 4 and Landsat 5 geolocated data was determined. Results of the correction effort 
show that the models significantly improved the Landsat 4 data and the initial result of 
the cross calibration effort looks very promising. 

A relative calibration was performed first using the IC. A data set of 260 scenes 
was studied. Due to a gap in the data, two models for each reflective band were required. 
Only the four reflective bands were studied at this point due to lack of an outgassing 
correction. The data was further restricted to only reverse scans, to stay consistent with 
Landsat 5’s IC study, and data from the primary lampstates [001], [010], and [100]. 

Once the data were selected and processed, the lamp data was recorded. This record 
included readings from each detector, in every band in the PFP, for lampstates [001], 
[010], and [100]. 

The goal of the IC lifetime analysis was to determine a single model for each 
band in the PFP. However, to arrive at this model, the lifetime trends had to be analyzed 
at multiple lower levels. First, detector trends were studied to verify no single detector 
was performing differently than the rest. Then ‘odd/even’ trends were studied. An 
odd/even trend refers to averaging all odd numbered detectors together and all even 
detectors together. This is a valid study due to the physical layout of the detectors within 
a band on the TM instrument. Next, the lampstate trends were studied. This trend is the 
average of all the detectors within a band, after any necessary corrections were made for 
odd/even effects, for each of the three lampstates. It was determined that lamp [001] was 
performing differently than the other two lamps; therefore, it was not included into the 
final model. Finally, all the data from lampstates [010] and [100] were averaged into 
band trends and models are determined. 

Now the band trends could be studied and lifetime IC models could be 
determined. It was decided to use two models for each band, one for the pre-gap data and 
one for the post-gap data because in each band there is a significant difference between 
the trends after the gap. These models are in the form of a Linear + Exponential 
regression. This model was chosen because of the work done previously on the Landsat 
5 TM. The model coefficients are displayed in Table 1 and the Goodness-of-Fit statistics 
are shown in Table 2. The RMSE values for each model are very good, less than 0.01 
RMSE. The R-squared values are also acceptable. Bands 1 and 4 are quite good, R 2 > 

0.7, but Bands 2 and 3 are a bit lower, around 0.5. 
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Table 1 : IC lifetime models for Landsat 4 PFP bands y = ae (bx) + cx + d 


Coefficients for IC models 
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Table 2: Goodness-of-Fit statistics for the IC lifetime models 
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After the IC models were created, the RelGain striping in the image was analyzed, 
modeled and removed. This effort was largely simplified by referring to the work 
previously performed on Landsat 5. This method produced a model that considered data 
spanning the lifetime of the instrument and produced improvements greater than 75%, 
assuming the scene being improved had some striping already present in the image. If 
there was no striping present (Striping Ratio before correction < 4) and the scene mean 
was extremely low (mean < 20), the models failed. In this case, the old method of 
histogram equalization correction should still be used. However, the number of scenes 
with both of these restrictions are very few; most of the scenes were improved drastically, 
many by > 90%. Figure 1 shows an example of RelGain and the results after correction. 



Figure 1: Relative Gain qualitative analysis of Peru B4; Day 267, 


Year 1988, Path 007, Row 064 Correction of 96.7%. 
A) Before Correction. B) After Correction. 






















































The last step of this effort was to attempt a radiometric cross calibration 
(CrossCal) was for the Landsat 4 and 5 TMs. To perform this study, Landsat 4-Landsat 5 
scene pairs were obtained and compared. An assumption was made that the ground and 
atmosphere, and therefore the radiance, did not change during the time between samples. 
These images were to be homogeneous and near perfectly geolocated. Hence, any 
difference in the image would be due to the difference in the instruments. This difference 
could be modeled and then removed from the data. To model these differences, eight 
small sub-images were chosen throughout each scene-pair studied; unfortunately, only 
two scene-pairs were obtainable at the time of the study. 

Once these pairs had been chosen, data points were calculated by determining the 
average DN value of a sub-image. The data was analyzed by plotting L5 data on the x- 
axis and L4 data on the y-axis and a regression line through these points, anchored at the 
origin, was determined. Also included on these plots are one-to-one lines for reference; if 
the reference line fell directly on the one-to-one line, no correction would be needed. 

An example of this is shown in Figure 2 and the gain results are displayed in Table 3. 

Also in Table 3 are the 1 -Sigma values of the residuals and the corresponding percent 
error. All bands had less than 2% error except band 5. 

Table 3: L4/L5 Gain Ratio results for each Band 


Band 

Gain 

iqsissbs 

Percent Error 

Band 1 

1.058 

2.3 

1.7% 


1.052 

1.3 

1.8% 

Band 3 

1.035 

1.6 

1.6% 


0.978 

1.5 

1.4% 

Band 5 

0.889 

6.1 

3.3% 

Band 7 

1.003 

1.9 

1.7% 
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Figure 2: L4/L5 Gain Ratio. Band 1 
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Band 5 was studied further and it was determined that much of the excess error 
was due to the outgassing cycle. While, at the time of this study, there was no outgassing 
model for band 5, the appropriate outgassing cycle had been plotted. By normalizing one 
data point to the other this error was reduced to 1 .7%, which was equal to the average 
error across the bands. 

These models produced CrossCal points on the specific days of the scene pairs. 
However, a lifetime model was desired. To do this, the IC models were scaled to the 
CrossCal points, resulting plots are shown below. The blue line is the IC model, the ‘x’ 
are the CrossCal points, and the red line is the adjusted calibration model. 


L4 absolute gain model, Band 1 14 absolute 9™ model, Band 2 
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Lastly models for the absolute gain of Landsat 4 were developed. The equation for the 
model is the standard exponential plus linear equation that was used with Landsat 5 TM: 

Gl 4 = Ae Bt + Ct + D 



Coefficents for L4 Absolute Gain Model 

Model 

A 

B 

c 

D 

pre-gap Band 1 

0.1565 

0.0210 

-4.565E-05 

1.428 

pre-gap Band 2 

0 

0 

1.913E-05 

0.679 

pre-gap Band 3 

0.3044 

0.0270 

1.357E-05 

0.912 

pre-gap Band 4 

0.0815 

0.0120 

2.608E-05 

1.021 

post-gap Band 1 

-0.1599 

-0.0004 

1.820E-04 

1.327 

post-gap Band 2 

-0.0001 

-0.0016 

2.888E-05 

0.643 

post-gap Band 3 

-3.903E-05 

-0.0019 

3.500E-05 

0.878 

post-gap Band 4 

-1.971E-05 

-0.0020 

3.181E-05 

1.004 


Summary 

A lifetime radiometric model for the absolute gain of Landsat 4 has been developed. The 
key features of this model are that it is based on use of the IC an interpolator over the 
lifetime of the instrument. The absolute gain of the instrument was obtained by cross- 
calibration with Landsat 5 using nearly concurrent scene pairs of invariant sites. A 
relative gain model was developed for the instrument, as well, to remove striping and 
radiometric errors due to non-uniform detector gains within a band. 
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Landsat 5 Scene Invariant Analysis 


Introduction 

During the 2003 contract year, work was begun on alternative methods to determine/ 
validate the absolute radiometric gain of the L5 TM sensor (bands 1...5 and 7) via 
examination of “invariant” sites. Initial results and methodology were reported in the 
2003 Annual Report based on processing of Yuma and RRV data sets. The early results 
indicated a measurable change in sensor response for bands 1 & 2 of approximately -12% 
and -10% respectively over the initial 6000 days of operation. Effort was expended to 
incorporate a linear decay term into the accepted IC based gain model to account for 
effects observed in the scene invariant data. The results were reported at the Calibration 
meeting hosted at SDSU. During 2004, additional work was undertaken to further process 
the Yuma data set and enhance the gain model for bands 1 and 2. Results were reported at 
the Lake Tahoe calibration meeting, and are summarized in this final report. 

Results 

Efforts were limited to examination of TM bands 1 & 2 where a reasonable amount of 
sensor responsivity change was observed above the noise level. The primary gain model 
extension attempted was to fit normalized scene observations with the following 
empirical model: 

^ ^ b + [c * exp (d * ( DSL - s))] 

This model allowed for initial and final state existence, along with an exponential 
transition between the states. Prior to fitting to the Yuma data set, a Savitzky-Golay filter 
was applied to smooth the data and account for non-periodic sampling. Results of the 
SDSU processing effort are presented in Figures 1 and 2. 
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While unexplained oscillations are apparent in both Band 1 & 2 data sets, it was decided 
to denormalize the resulting model, tie the results to the 1999 L7 cross-calibration and 
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compare to the accepted IC and vicarious gain results. Figures 3 and 4 present the final 
models that were developed. Shown are two potential scene invariant gain models which 
differ due to the input data utilized (The curves labeled SDSU Invarient was based on 
large scene regions, while the curves labeled GSFC Invarient were modeled with small 
region means provided by Ed Kaita at GSFC). 
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Examination of the resulting gain models, in comparison to the IC gain model and 
vicarious data points, reveals some unexpected results. First and foremost, a significant 
difference exists in the models wherein the accepted exponential responsivity decay 
during the early days of the mission is not evident in the Yuma scenes. Secondly, a 
transition was observed in the Yuma data approximately 3000 days after launch, but not 
observed in the IC data. The vicarious data does little to resolve the discrepancy, splitting 
the difference for Band 1 IC and scene invariant models, and indicating a higher gain 
than either the IC or scene-based models for Band 2. Critical observations during the 
1989-92 period which might help to resolve the matter are not available. The end 
conclusion is that insufficient evidence was provided by the scene invariant approach to 
justify changing the status quo regarding how DN’s are converted to radiance for Level 1 
L5TM products. Development of the TMIAS by NCEROS provides an opportunity for 
further study of both the IC and scene based methods. 
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SUMMARY 


The previous sections have outlined work done in six major categories. 

First, the Landsat 5 IC has been extensively studied and trended. The major change that 
has been observed in the past three years is a significant degradation in the performance 
of all lamps. It is unlikely that the lamps will be useful after about 2002 for any 
radiometric trending. Fortunately, the response of the instrument has been stable for 10 
years prior to this occurrence suggesting that, although the lamps are failing, there is 
reason to believe the detectors will remain stable. 

The second major area of study is vicarious calibration. The past three years has seen a 
significant change in personnel and capabilities. After a few years of not being able to 
perform a full vicarious calibration, the SDSU group is again able to provide absolute 
gain estimates using ground-based data collection and atmospheric modeling. In 
addition, an emphasis has been placed on accurate, repeatable data collection, as well as 
being able to provide calibration estimates in a ‘production environment’ with a rapid 
response time. Gain estimates have been comparable with results obtained by other 
investigation teams. 

One of the more troubling aspects of Landsat radiometry has been with respect to the 
relative gains of detectors within each band. This has historically resulted in a striping 
pattern on imagery that is not only annoying but a source of radiometric error. Histogram 
equalization using data from individual scenes has been the heritage method of 
addressing this problem. A study was undertaken to expand that method using multiple 
scenes over the lifetime of the instrument to develop a more accurate model of relative 
gains. This approach relies on the stability of the detectors relative one another as has 
been borne out on many image examples. Using the results of this study, relative gain 
correction is as simply as looking up the relative gain of each detector and applying a 
simple multiplicative factor for correction. 

In the cold focal plane, the outgassing effect has continuously plagued radiometric 
models for Bands 5 and 7. Work under this project has developed a Landsat 5 outgassing 
model that treats the effect as a thin film phenomenon. Based on this model an accurate 
and robust correction algorithm has been developed that essentially removes this 
radiometric error. This technique has been implemented into standard processing at the 
EROS Data Center. 

An absolute gain estimate for Landsat 4 has been developed. It is based on the data 
collected by two lamps of the IC. These data sets are tied to a cross-calibration with 
Landsat 5 using a pair of scenes collected at nearly the same time and over the same 
region. Results from this study indicates a cross-calibration with L5 that is accurate to 
around 3%. Unfortunately, gaps exist in the L4 data sets and it is difficult to extend this 
model accurately back to the day of launch. However, it is felt that to date this is the 
most accurate absolute gain estimate of Landsat 4 available. 
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The last area of study involved using invariant scenes as an aid in understanding the 
absolute calibration of Landsat 5. This approach indicated very consistent results with 
Bands 3-5, and 7. However, Bands 1 and 2 results were different during the earlier 
portion of the Landsat 5 lifetime. Unfortunately, during the critical part of the lifetime 
analysis, where a large change is indicated, there is very limited data available to study 
the effect closely. Thus, initial thoughts are to disregard this approach until more data is 
available and a more thorough study can be performed. 

Lastly, it can clearly be seen that several major contributions to the understanding of the 
radiometric performance of Landsat 5 have been achieved by this group. This is 
especially important in light of the recent failure of L7 and the need to continue data into 
the future. 
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